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ABSTRACT 

This  paper  considers  an  alternating  renewal  process 
with  time  to  failure  being  distributed  exponential  (A) , 
and  with  a  general  distribution  for  time  to  repair. 
Bounds  and  the  precision  of  the  bounds  for  the  avail- 
ability function  were  obtained. 

A  computer  program  was  written  to  solve  for  the 
availability  function  and  some  other  quantities  when  the 
down  distribution  was  a  gamma  (a, 3),  a  integer.   The 
availability  function  with  down  distribution  being  a  mix- 
ture of  gamma- functions  is  considered. 
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I.   DESCRIPTION  OF  ALTERNATING  RENEWAL  PROCESS 

A.  SYSTEM 

The  system  considered  consists  of  a  single  unit  or 
component  which  is  repaired  upon  failure  and  then  re- 
turned to  operation. 

In  the  model,  the  following  assumptions  are  made: 

-  Time  to  failure  is  a  random  variable  (R.V.)  X  with 
distribution  F  (specified  below) . 

-  Time  to  perform  a  repair  is  a  R.V.  Y  with  dis- 
tribution G  (specified  below) . 

-  Repair  commences  immediately  upon  failure. 

-  Unit  is  returned  to  operating  state  immediately 
upon  completion  of  repair. 

-  Unit  is  a_s  good  as  new  after  repair. 

-  Consecutive  operating  times  between  failures  (as 
well  as  repair  times)  are  independent. 

Thus  the  system  can  be  modeled  as  an  alternating  re- 
newal  process . 

B.  DISTRIBUTIONS 

The  time  to  failure  -  termed  uptime  -  is  throughout 
assumed  to  be  exponentially  distributed,  with  probability 
density  function  (p.d.f.): 

(1.1) 

f(t)  = 


Three  different  distributions  for  time  to  repair 
termed  downtime  -  were  treated: 

-  Gamma  distribution. 

-  Mixture  of  two  exponential  distributions. 

-  Mixture  of  gamma  and  exponential  distribution. 
The  respective  p.d.f.'s  are: 
(1.2) 

0  ;   t  <  0 

3  >  0 

a  =  1,2,3, • • • 


t*-1  e"et  ;   t  >  0 


g(t)  = 


(1.3) 

g(t)  = 


m31e"eit  +  (l-m)32e"e2t    ;   t  >  0 
0  ;   t  <  0 


BlfB2  >  0 
0  <  m  <  1. 


t  >  0 
t  <  0 


C.    LAPLACE  TRANSFORMS;  GENERAL 

Denoting  the  operating  state  by  1  and  the  failed  state 
by  0,  the  objective  is  to  find: 

-  Plx(t):   The  "availability"  of  the  system,  i.e., 
the  probability  that  the  system  will  be  "up"  or  operable 


when  called  upon  at  a  random  time  t,  given  that  it  starts 
in  state  1  at  t  =  0 . 

The  following  are  also  of  interest: 

-  Mxl(t):   the  expected  number  of  visits  to  state  1 
in  period  (0,t),  given  that  the  system  is  in  state  1  at 
t  =  0. 

-  P   (t) :   probability  system  is  in  state  1  at  time  t 
given  it  starts  in  state  0  at  t  =  0 . 

-  EjjCT):   expected  time  system  is  "on"  in  (0,T). 
It  is  shown  in  [1]  that 
(1.5) 


Mn(t)  =     M01(t-x)dF(x) 


(1.6) 

M 


o,(0  =  f  [1  +  Mix(t-X)]  dG(x) 


Taking  Laplace  transform  on  both  sides  of  (1.5)  and 

(1.6),  denoting  the  transform  by  circumflex  ("hat"),  i.e. 

M..(s)  =  P      e~st  M.-(t)dt,  one  obtains: 
1-)      o         1-1 

(1.7) 

f(s)-g(s) 


Mlx(s)  = 


s(l-f(s).gts)) 

In  similar  manner: 
(1.8) 

p1iW  -       l^isL, 

s(l-f(s).g(s)) 


(1.9) 

p  ro  -  Ci-f(s))  g(s) 

r 0 1 vbJ 

s(l-f(s)-g(s)) 


(1.10) 


E(s)  ■      X-f^ 


s2(l-£(s)g(s)) 


Only  for  some  of  the  distributions  considered  can  the 
inversion  be  accomplished  in  a  closed  form.   For  the  gen- 
eral case  with  gamma  down  distribution,  a  computer  program 
was  written  to  give  numerical  solutions. 

If  a  digital  computer  is  not  available,  good  approxi- 
mations are  of  interest.   These  -  together  with  bounds  - 
have  been  developed  for  the  availability  function. 

Some  solutions  in  closed  form  in  the  simpler  cases 
are  provided. 


II.   EXACT  SOLUTION  -  LAPLACE  TRANSFORM 


Laplace  transforms  of  the  up-and  downtime  probability 
density  functions  considered: 
(2.1) 


f<S)  =  T^ 


A+s 


(2.2) 


(s)  ■ 


;  +  s 


(2.3) 


(2.4) 


g(s)  =  m  g— - 


1-m 


+  s 


Ref.  (1.1) 


Ref.  (1.2) 


Ref.  (1.3) 


Ref.  (1.4) 


In  the  case  of  (2.2)  one  obtains  upon  substituting  into 
(1.7)  to  (1.10)  and  rearranging: 
(2.5) 


MM(s)  = 


X6V 


s-q(s) 


where , 


(2.6) 


q(s)  =  (s  +  3)a  +  X  Z   B^^s  +  e)01"11   . 
i  =  l 


(2.7) 


Pn(s)  = 


(s  +  3) 
s-q(s) 


(2.8) 


Poi(s)  = 


q(s) 


(2.9) 


"i»C»)  " 


(s  +  B) 


In  cases  (2.3)  and  (2.4),  only  Plx(s)  has  been  con- 
sidered: 

(2.10) 

(s+Bx)  (s  +  62) 


PM(s)    = 


(2.11) 


PuCs)    ■ 


»[sa+(X+fi1+0a)s+eiCA+Ba)-^aCB1-e   )] 


(s+sj   (s+e2) 


s[s3+As2+Bs+C] 


Ref.    (2.3) 


;    Ref.    (2.4) 


where 


A  =    Bx    +    2B2    +    A 

B   =    B2    +    2B2(A+BJ    +    XBxin 

C   =    XB2(l-m)    +    B^2    +    2AB    B  ra 


III.   INVERSION 

A.  COMMENTS 

Consider  the  transform  of  a  rational  function  P(s)/sq(s) , 
where  the  degree  of  P(s)  is  less  than  sq(s) .   This  is  a  suf- 
ficient condition  for  the  use  of  the  residue  theorem  in 
finding  the  inverse  Laplace  transform  [2],  thus: 

(3.1) 

F(t)  =  E  (Residues  of  est  f(s)  at  poles  of  F(s)). 

As  is  shown  in  Appendix  B  only  simple  poles  exist  if 
£  >  X.      This  simplifies  the  inversion. 

Equation  (2.9)  has  a  double  pole  at  s  =  0 .  The  ap- 
propriate adjustment  to  the  procedure  used  is  considered 
later. 

In  the  sequel,  primarily  the  inversion  of  Plx(s)  will 
be  considered.   However,  since  the  method  used  is  applicable 
to  all  cases  where  Laplace  transform  is  of  type  (3.1), 
the  discussion  in  some  detail  of  inversion  of  Pu(s)  will 
show  how  the  other  results  were  obtained. 

B.  INVERSION  FORMULA 
Referring  back  to  (2.7) 
(3.2) 

Pll(t)  =  L-l  (lilil"] 
\sq(s)  | 

q(s)  is  a  polynomial  of  degree  a.   It  is  shown  in  Appendix 

A  that  all  roots  are  complex,  except  when  a  is  odd,  in 

which  case  one  of  the  roots  is  real  and  negative. 
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where 


Factoring  sq(s)  into  linear  terms: 
(3.3). 

sq(s)  =  s(s-rj  (s-r2)  ••  •  (s-ra) 

r.  =  a.  +  b.i      (b .  =  0  if  a  is  odd) 
J     3     J         * 

r .  ,  =  a.  -  b  .i      j 
Using  the  residue  theorem: 


Pn(t)  =  E  Residues 
a+1 


/         lim 


(s-rk)(s+3)a   est 
sCs-rJ  (s-r2)  •••(s-rk)  •••  (s-ra) 


C3.4) 


aX 


■a      *\t 


k=l 


£  -   1,2, •••,<*      £  ^   k 


Denoting 


C.VB)      -   vk  ♦  wki 


n 

k=l 


rk(rk"rP  =  xk +  >v   l* k 


then 


(3.5) 


Plx(t) 


azr-T  +    E((vvx,  +wlr y,  )cos   b,  t   +    (v,y 
3  +  aX         ,  v     k   k      k'  k'  k  k* 
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(3.5)Cont'd. 


e3kt 

sin  b,  t) 

K       2.2 

k  =  1,3,5, ••• ,a-l. 

Incidentially ,  for  a  bounded  solution  to  exist  the 
real  part  of  all  roots  must  be  negative. 

The  numerical  solution  is  contingent  on  finding  the 
roots  of  sq(s).   This  can  be  done  algebraically  for  a  <  4. 
The  method  used  is  outlined  in  Appendix  C.   In  addition, 
for  a  =  2  the  inversion  of  (3.5)  in  explicit  form  is  shown  in 
Appendix  D  to  illustrate  the  method. 

For  a  polynomial  of  unrestricted  degree  the  computer 
program  described  in  Section  V  was  written. 

If  a  computer  is  not  available,  approximating  formula 
P^Ct)  are  given  in  Section  IV. 

For  E1X(T),  equation  (2.9),  the  double  root  at  s  =  0 
leads  to: 

E   (T)  -  lim  J-     s2(s  +  S)a  ^ 

s-0  ds   s2(s-rI)(s-r2)---(s-ra) 


I 


(s-r*)(s+B)aesT 
lim  5 


k=1   s+rk   a 

hi 

which  upon  derivation  and  letting  s  ■*  0  lead  to 


(3.6) 


E   (T)  =  -  ^Ce+Wcx-l))  +   X6   T 


(6+aX) 


aX 


I 

k=l 


lim 


(s-rv)(s+6)a  esT 


*k  sz(s-r1)(s-r2)«-'(s-rk)---(s-ra) 


A  similar  procedure  as  outlined  for  Plx(t)  is  used  for 
the  last  term  (modified  by  multiplication  by  s) . 
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IV.   APPROXIMATIONS  AND  BOUNDS 


A.    GENERAL  DOM  DISTRIBUTION 
Consider 


s(l-f(s).g(s)) 


where  uptime  X  -  exponential,  downtime  Y  ~  G,  unspecified 
Then, 

1--*- 


sjl-^-gCs) 


(4 

.1) 

^ 

Pn 

(s) 

d 

ivision: 

(4 

.2) 

[i-xntoi] 


Observing  that  l-g(s)/s2  is  the  Laplace  transform  of  the  equi 
librium  distribution  of  G  (Reference  4)  except  for  a  factor 
1/E[Y],  then 
(4.3) 


P11(s)  =  1  -  RGe(s)  +  R2(Ge(s))2  -  R3(Ge(s))3+... 
where 
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With  R  <  1,  convergence  is  ensured.   This  implies  that 
Pu(t)  >  \   as  t  ■*>  °°  (steady  state  result).   Inverting  (4.3) 

Pxl(t]  =  1  -  RGe(t)  ♦  R2G^2)(t)  -  R3G^3)(t)  +••• 

where   G^1J(t)  denotes  the  ith  fold  convolution  of  G  with 
itself. 

Letting  t  ■*  °°, 

limPllCt)  -  1  -  R«  R*  -  R»  ♦  •..  -  E[xfI*}Y] 

one  obtains  the  steady  state  result. 

B.    BOUNDS  AND  PRECISION 

It  is  a  reasonable  proposition  to  try  to  invert  the 
first  term  in  (4.3).   Inversion  of  the  other  terms  may  be 
difficult  unless  G   is  a  simple  distribution. 

Even  though  (4.3)  is  not  a  series  in  the  time  domain, 
by  probabilistic  considerations  th.e  effect  of  truncation 
may  be  deduced. 

Noting  that  R1  G^  (t)  >  R1  +  1  G^1  +  1)  (t)  ,  then  the  ap- 
proximation of  P1]t(t)  will  give  a  lower  or  upper  bound  on 
Pxi(t)  depending  on  whether  the  last  term  is  negative  or 
positive ,  e . g. :  » 

(4.4) 

Pxl(t)  >  1  -  R'Ge(t) 

(4.5) 

P„(t)  <  1  "  R'Ge(t)  +  R2Ge(2)(t). 
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The  precision  with  which  the  bounds  are  established 
can  be  obtained. 

Suppose  one  truncates  after  one  term.    Then  the 
maximum  value  of  the  remainder  is  less  than  R2 .   Thus, 

(4.6) 

1  -  RGe(t)  +  R2  >  PM(t)  >  1  -  RGe(t). 

Having  found  the  lower  bound  for  Plx(t),  then  this  value 
is  accurate  to  within  R2 .   In  a  rather  typical  case  with 
R  =  1/10  (i.e.,  steady  state  =  0.909),  the  lower  bound 
is  established  to  within  one  per  cent. 

If  it  is  feasible  to  compute  the  upper  bound  by 
taking  two  terms  in  the  expansion,  then  -  by  similar 
argument  -  the  precision  of  the  upper  bound  will  be  R3 , 
i.e.,  to  within  1/10  per  cent  in  the  case  considered.   With 
R  relatively  large,  the  precision  is  poor.   If  feasible, 
more  terms  must  be  computed  for  a  desired  precision. 

C.    GAMMA  DOWN -DISTRIBUTION 

Consider  case  with  Y  *  Gamma  (a, 3)  a  =  1,2,3, •••: 
(4.7) 


a- A    1  - 
R-G.(s)  - 


(3+s)a-ea 


=  X- 


3      s2(£]        s2(3+s)a 


Noting  that  s  =  0  is  root,  one  may  simplify  (4.7)  as  fol- 
lows : 

(3  +  s)a-3a  =  3^3a-1s^a"2(3^s)s  +  3a"3(3  +  s)2s+>-»  +  (3  +  s)a-1s-ea 
s2(3+s)a  s2(3+s)a 
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1  1 

s  e 


6+s  J  [3  +  sj 


a-1 


B  +  s 


(4.8) 


a 


Il^l1    • 


RVS)  -H   y 


Substituting  into  (4.3) 
(4.9) 


Pxl(s)  = 


1 


Each  term  within  the  bracket  contains  gamma  trans- 
forms, making  the  inversion  simple. 

As  an  example,  consider  case  with  a  =  2.   Then  (4.9) 
becomes : 

(4.10) 


x      b  +   b 
b  [b+s   \e+s 

B  +  s      B  +  s 


X 


+  s 


s-PM(s)  ■  1 


+  2 


The  resulting  inverse  for  first  approximation: 
(4.11) 

Pn(t)  >  1  -  |  (1  -  e"8t  +  1  -  e"6t  CBt+1)). 

Using  the  correspondence  between  poisson  and  gamma  dis 
tributions,  the  value  can  be  obtained  from  cumulative  pois- 
son tables. 


17 


D.    MIXTURE  OF  GAMMA  DOWN- DISTRIBUTION 

Consider  a  mixture  of  two  gamma-distributions  with 
parameters  (K,3),  (&,y)  and  weighted  (1-m)  and  m  respec- 
tively.  The  first  term  in  the  expansion  becomes: 

(4.12) 


PnOQ  - 


1 


1 


i^Hl(s)  +  *H2(s) 


where 


HiW  =  3S7  +  ^ 


J3_ 


H,(s) 


-X-  +  LL 

Y  +  s    y+s 


The  first  term  is  seen  to  be  the  linear  combination  of  the 
two  terms  in  the  down-distribution.   For  case  K  =  2, 
£  =  3,  the  inversion  is: 
(4.13) 


Pn(t)  >  1 


where 


~^     Mt)  +  J  H2(t) 


Hx(t)  =  1  -  e"3t  +  1 


e"et  (St+1) 


H2(t)  =  1 


Yt 


1  -  e 


Yt 


(Yt  +  D  +  1 


-yt 


^  +  Yt  + 


In  case   more  terms   are   required,   the   inversion 
will  be  somewhat  complicated  in  that  there  will  appear 
cross-terms  between  the  two  distributions.   The  precision 
of  (4.13)  will  be 

' (1-m)2  +  mf 
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E.    WEIBULL-DISTRIBUTION;  SPECIAL  CASE 

A  lower  bound  can  be  obtained  when  downtime  is  dis 
tributed  Weibull  (3,2).   The  reliability  function  is: 

(4.14) 

G(t)  =  e"6t2  ,  3  >  0. 
The  expected  value  is: 

bit]  =  vy  . 

The  equilibrium  distribution: 


dx 


Ge(t)  = 


Em 


which  becomes  upon  re-arranging: 


-  r<o  _VT 


t  "I 

f  ,"^x2(23)  ,   I 


V^V^  j 


1 

■  -fa 


(4.15) 


v«  =^1 J  ^: 


(23) 


dx. 


Thus,  G  (t)  is  -  normal  (0,1/23).   The  lower  bound  becomes 
(4.16) 

Pll(t)  >  1  -  xVf 


where  Z(t//2~3)  is  normal  (0,1).   Thus  standard  normal 
tables  can  be  used  for  evaluation.    The  corresponding 
precision  is  A2*tt/3. 
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1 .   Weibull -Distribution;  Decreasing  Failure  Rate 

For  the  special  case  that  the  Weibull  distribution 
has  shape  -  parameter  K  =  1/p,  p  =  1,2,3, •••,  the  lower 
bound  with  one  term  in  expansion  can  be  found.   The  equi- 
librium distribution: 

(4.17) 

t 
e 
o 


G  (t)  -  4   f 

e^  J        p!   J 


1 

zp  dz 


which  upon  substitution  becomes 


(4.18) 


Ge(t) 


.r_£_ 

J        Cp-DI 


-ex    D-i 


xF    *   dx 


where 


Equation  (4.18)  is  the  gamma-distribution  in  x(=zp)  with 
parameters  (p,3).   The  lower  bound  becomes: 


Pn(t)  >  i  -  x   ^  r  (t;p,3) 


with  precision 


El  A 

oP 


As    an   example,    let   k   =   h,    i.e.,    p    =    2.      Thus 


P»Ct]     >    1    -    ^    at 


1    -    e 


/t 


/t   +   1 


)]• 


The  value  of  term  in  brackets  can  be  obtained  from  poisson 
t  ab  1  e  s  . 
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V.   COMPUTER  PROGRAM 

The  computer  program  was  written  in  FORTRAN  IV  and 
run  on  the  IBM  system/360  model  6  7  computer  at  the  U.  S. 
Naval  Postgraduate  School.   The  program  listing  is  shown 
in  Appendix  E. 

The  main  sequence  in  the  program  is: 

(i)  Input  data:  a,B,A,  time  increments  and  number 
of  time  steps  (named  respectively  M,  BETA,  UPRATE,  TSTEP, 
L  in  program) . 

(ii)   Generate  coefficients  of  polynomial  in  denom- 
inator of  (2.7). 

(iii)  Call  library  subroutine  DPOLRT  to  solve  for 
roots,  using  Ne w t on - Rap h son's  method. 

(iv)   Compute  and  store  coefficients  for  cosine, 
sine  terms  in  Pn(t),  P01(t),  Mu(t)  and  Elx(t);  e.g., 
see  (3.5)  and  (3.6) 

(v)    Sum  cosine  coefficients  for  Pxl(t)  and  add  to 
steady  state  value  of  Plx(t),  to  give  an  indication  of  ac 
curacy  of  final  result.   Theoretically  the  sum  should  add 
to  1.   Flagging  value  IR  provided  by  subroutine  to  in- 
dicate accuracy  of  root  extraction  should  read  zero. 

(vi)   Enter  DO-loop  for  incrementing  time  and  for 
final  calculations. 

(vii)  Print  results:    t,  Pn(t),  P0i(t),  Mn(t), 
Eu(t),  named:   T,  PROBUP,  DPRBUP ,  RNWAL ,  EXPUP  respec- 
tively in  the  program. 

Sample  output  is  shown  in  Appendix  F. 
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APPENDIX  A:   NATURE  OF  ROOTS 

Polynomial  considered: 
(A.l) 

P(x)  =  xa  +  XCx01"1  +  3xa"2+...+  3a_1)  ;  B  >  A  . 

Multiply  both  sides  by  x-B,  to  obtain: 

(A. 2) 

Q(x)   =  xa+1  -  (B-X)xa  -XBa  =  P(x)(x-B). 

Let  a  be  even. 

Then  by  Descartes'  Rule  of  Sign  [3],  with  only  one 
variation  of  sign,  there  exist  at  most  one  real,  positive 
root  (clearly  x  =  B) . 

Consider 

(A. 3) 

Q(-x)  =  -|x|a+1  "  (B-X)xa  -  XBa. 

Since  there  are  no  variations  in  sign,  no  negative 
real  roots  exist.   CONCLUSION:   For  a  even,  (A.l)  can  only 
have  complex  roots. 

Let  a  be  odd. 

By  (A. 2),  one  variation  indicates  the  positive  root 
X-B.   Consider 

(A. 4) 

Q(-x)  =  xa+1  +  (B-A)|xa|  -  XBa. 

One  variation  imply  one  real,  negative  root.   CONCLUSION: 
For  a  odd,  (A.l)  has  one  negative,  real  root,  (a-1)  complex 
roots  . 
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The  roots  of  (A.l)  are  of  the  form 

r."  =  a.  +  b.i   ;  j  =  1,2, •••a  with  bt  =  0  if  a  is  odd 
Since  x  =  s+3,  the  roots  of  s  will  be 


r!  =  r,  -  3  -  a .  -  6  +  b.i 


and 


£  r!  =  -  X  -  a3. 

j  =  l 

That  all  a.  <  0  can  be  verified  by  Hurwitz  Criterion  (3) . 
However,  since  E[X+Y]  <  °°  and  F*G(t)  is  non-lattice,  the 
system  will  be  stable,  implying  that  a.  <  0  for  all  j  [4] 
(*  denotes  convolution) . 
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APPENDIX  B:   NON-EXISTENCE  OF  MULTIPLE  ROOTS 

Polynomial  considered: 
(B.l) 

P(x)  =  xa  +  A  (xa-1  +  3xa"2  +---  +  |3a"1)  ;  3  >  A. 

For  simplicity,  let  a  =  6,  and  assume  that  a  double  and  a 
single  complex  root  exist: 

rl,2  =  a  ±  bi 

r3,4  =  a  ±  bi 

r5,6  "  *  ±  di. 
Consider  the  usual  relationship  between  roots  and  coeffi- 
cients [3]  : 
(B.2) 

4a  +  2c  =  -X  ;  sum  one  root  at  a 

time 
(B.3) 

(a2+b2)2(c2+d2)  =  A35  ;  sum  6  at  a  time 

(i.e. ,  product  of 
(B.4)  roots) 

2(a2+b2)+ (c2+d2)+4a2+8ac  =  3A   ;  sum  roots  two  at  a 

time . 

Using  Cauchy's  inequality  on  (B.3): 

(B.5) 


2(a' 

"(a'+b^'Cc'+d*) 


+b2)+(c2+d2)  >  *f 
By  (B.4) ,  then 
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(B.6) 

2(a2+b2)+c2+d2  =  BA  -  (4a2+8ac) 

=  BX  -  4((a+c)2-c2)  <  BA 

since  (a+c)2-c2  >  0. 

Equations  (B.5)  and  (B.6)  are  in  contradiction,  since 
3/X35  >  3X  for  all  values  of  3  >  X.   Thus  double  roots  are 
not  possible.   The  proof  can  be  extended  to  any  a  and  any 
multiplicity  of  complex  roots  considered.   The  proof  holds 
also  for  a  odd. 

The  two  inequalities  established  for  a  =  6  hold  in 
general  since  all  terms  of  form  (a2+b2)  will  appear  in 
isolation  from  all  possible  cross  terms  of  form  (a»c). 
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APPENDIX  C:   POLYNOMIAL  DEGREE  3  AND  4 


Provided  3  >  A,  the  roots  of  polynomial  of  degree  3 
and  4  can  easily  be  found  [5]. 
a  =  3 

P(y)  =  y3  +  Ay2  +■  ASy  +  A32  ;  y  =  s  +  3. 

Denote : 

a  =  |  (3A3-X2) 

b  =  27  (2A3-9A23+27A32) 


A  =    -  %   +  Vb2  A  a3 


4    27 

M 


-*-^! 


Then  the  roots,  solving  for  s,  are: 

r1=A+B-3  -  t?  (real  root) 

r2  3   =  -  |  (A+B)  ±  |  (A-B)  YTi-3-^  (complex) 

a  =  4 

P(x)  =  x-  +  Ax3  +  A3x2  +  A32x  +  A33  ;  x  =  s 
Denote : 


a  =  |  (A-63)X32 
b  =  ^|  (93-4A)A233 
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—  i-»r 

a.] 

a»m| 


38  =  arccos 
y  =  m  cos  6  +  AS/3 


I^t  -  X3  +  y 


D1  2  =  |  I  X2  -  Xg  -  y  ±  (4X2B-8X32-X3)/4R|^. 

Then,  the  four  complex  roots,  solving  for  s,  are: 


rl,2  ■  -  I  R  "  (tX  +  S)  *  I  V 
r3>4  =  \   R  -  (Ix+6|  *  |  Dli 


The  simplicity  in  form  of  polynomial 
P(x)  =  xa+1  -  (g-X)xa  -  X£a 


make  solving  for  roots  for  a  >  4  by  manual  methods  easy 
for  example,  by  Graeffe's  Method  [3]. 
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APPENDIX  D:   AVAILABILITY  FUNCTION  IN  EXPLICIT  FORM;  a  =  2 

In  order  to  illustrate  the  inversion  method  used,  con- 
sider P„(t)  for  a  =  2: 
(D.l) 


Pxl(t)  =  L 


where 


(s  +  S) 


;((s  +  S)2  +  A(s  +  B)+X3) 


-  L 


-1 


(s  +  B) 


s(s-rj) (s-r2) 


t    =  \   (  -X-2B±/4AB-X2  i)   ;  B  >  A 


=  |  (-X-2B±Ai) 


;  A  =  /4XB-A2   . 


Then, 


sfs  +  812  est    ,.   fs-rj  (s+6l 2  est 

,  ,  CO  =  lim  — — - — ££- t  +  ±im  '  rrk2  >  ^  >~  -■.  ^ 

s(s-r.)  (s-r2)        sCs-r^Cs-rjj 


s-0 


+  lim 
s-*r, 


(s-r2)(s+B)2  est 
sCs-rJ  (s-r2J 


B2      (r^B) 


r,t 


(r2+B) 


+  7 r  e  l       +  7 =f 

rxr2   r1(r1-r2)        r2(r2-r1) 


r,t 


Since  the  last  two  terms  are  complex  conjugate  numbers, 
they  may  be  collected  and  expressed  in  terms  of  trigono- 
metric functions: 

B2 
PM(t) 


+  e 


i+|)t 


B2+2XB 

At  --t 

(r1^B)2r2(r2-r1)e2  *    (r2+B) 2r , (r , -r2) e 

r1r2Cr1-r2JCr2-r1j 
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8  +  2A 


+  2e 


-(B+j)t   (ac+db)cos^t  +  (ad-bc)sin^t 
c2+d2 


where 

C^  +  8)2  =  a  +  bi 

TliTl-T2)     =     C  +  di 

Substituting  in  actual  values,  the  result  is: 

8       2A 


(r2+6)2  =  a  -  bi 


r2Cr2-  rx)  =  c  -  di . 


pnW  =  *T: 


\A   cos  I 


t  +  (8 -A)  sin  yt  f  e 


C3+J)t 


2A   A(3+2A) 
Closed  form  solutions  were  also  found  for 


-(84)t 

=  _  A(28  +  A)  +  _A£_  t  +  Ae J  (282+ A2) sin  yt 

(B+2A)2    8+2A       (6+2A)2A 


A-  (28  +  A)cos 


?.,W 


H- 


8+2A    (8+2A)A 


(8+T)t 


j  (28  +  A)sin  |  +  A  cos  |t  5-  e 
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APPENDIX  E:   COMPUTER  LISTING 


I^Pi  TCTT    sc4L*8(*-H,n-Y) t0 OMPLEX*l 6 ( Z) 
DIMENSION    xrpFdOO  ),C  OF  dOC),  ROOT  RdOO),  ROOT  H  100) 
DTMENST  ON    ZR<  IOC)  ,  ZROOTdOO)  , ZNUM(IOO),  ZDNMdOO  ) 
DIMENSION    SCOEF(  100).GC0EF<  100) 
DIMENSION    Q,=  F=(100) 

DIMENSION  OOCOFFdOO)  ,0SC0FFU00)  ,OQcFcd00) 
DIMENSION  ROEFF( 100  ) ,  P.QCOEF  (  100),RSCPEF< IOC) 
DIMENSION  ^CrO^FClOO)  ,  ESCOEF(IOO)  ,50cFf=(100) 
DIMENSION    XMUM( 2)» XDMM( 2) ,XRDNM(2) tX^ (2) 

EQUIVALENCE    (ZZDNM,XpNM)t ( ZZNUM, XNUM) , ( ZR ONM, XD DNM ) , ( ZZR , XP ) 
M  =  5 

UPRATE=1.0nc/4.5D0 
R=TA=1C.0D0 
TSTEc=C.Cl°r 
C  GENERATE    POLYNOMIAL    AND    SOLVE    FOR    ROOTS 

J  =  M+1 

XCOF(J)=1.000 

XCOF(M)=UPRATF 

suvi=o.  cro 

03    1000    M=2.M 

IUP*IUP+1 

L=M-M+1 

XCOF(L)=UPPATE*BETA**IUP 
1000    CONTINUE 
r  '    mov    en?     TMVC°SIPN    HP    LAPLACE    TRANSFORM 

CALL    DPCLRT(XCOp»C0F,Mf ROQTR »POOTI»iFR ) 

DO    1500    T-l," 

ZPOPT(  T  )=nc»'FLX(DnOTP(  I  ).POOTId)  ) 

Z0(T)=7Ff-r-T(I  J-BETA 

ZMUM(I »=ZROOT< I )**M 
15C~>    r^NTTMl|£ 
C  ' "ABivc    DID    CALCULATE    NUMERATORS.       NOW    FOR    DENUMERATORS 

pn    2000    t=1,m 

»«=(I  ,M':.l  )    GO    T^    2001 

ZDO0P=ZR(l)-ZP(2  ) 

IFCM.GT.2)    r-P    T°    2002 

ZDNMd  )=ZPnrn*ZR  (1  ) 

£<■)    to    ^5^2 

2001  ZPROC=ZP(I)-ZR( 1 J 
I PR0D=2 

IF(M.F0.2)     GO    TO  -2501 
GO    Tn    2003 

2002  nYM'1Y=0.0 
IPR0P=3 

2003  PUM=0.0 

DO    2  5 OC     J=T  PROP»M 
IF((J.F0,I).ANd!(J.EO.M))    GO    TQ   2501 
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Ic<  J. EC!)     GO    TO    2500 
ZopnD=ZPPOO*(ZR< I)-ZM J)) 

2  5  00  CCNTTNUC 

2501  Z0NM(I)=ZPRGD*ZR(I1 
20 CO  CONTINUE 
C  K'~XT  P^llqws  SINE  AND  COSINE  COEFFICIENTS 

!CHFCK=CV2)*2-M 

IPdCHFCK)     3501,3502,3501 
3502    IM&G=I 

GO    Tn    3503 
3501     I«AG=2 

3  503    DUMMY*].. 0 

ILTm=m-i 

no   350^    !=IMAG,ILIMf2 

Z7PN^=Zr^M(  T  ) 

ZZNUM=ZNUM(  I  ) 

ZCQN»  =  Z^  M(  T>*  1°  (I  ) 

0MUM=x^U^( 1)*XPNM( 1)+XNUM(2)*X0NM<2) 

0DNM=XPNM( 1)**2+XDN*< 2)**2 

0r  OF-  (  T  )  =2  •  0*^0*  omijm/odNM 

SNUM=XfMJM(  1)*XDNM<  2)-XNUM(  2)*XPNM<1) 

SCnEc(  t  )=2  ,0D0*SNUM/ODNM 

DnrOFF(  !  )  =  2  .00  3*Pr"T^**v*XDMM(  1)  /ODNM 

CSfOEF(  I  )=2.0DO*«5TA**M*XDNM<  2)/0DNM 

RDNM=XRDN'M(1  ) **2+X RDMM (2)**2 

ROCPF  =  (  I  )  =2.  0?  0*UPR £TE*FET?"*  M*x0PNM(  1  ) /P^NM 

p^rn'rc  *  f  }  =  2    OOO^'JPR  AT  c*D'-~r  -  ~^  ~^M  ^X'-'  ""^I M  {  2)  /s  n  mm 

PONUM=:XNMM{l)*XRPNM(l  >+XNU:M(2>*XRDNM(  2 ) 

ESNUM=XNIJM(  1>*XRDNM(  ?)-XNUM(  2)*XRDNM(1) 

PSCPEF( I 1=2«000*PSNUM/RDNM 

EOCPSFt T  >=2 .ODO*EONUM/RDNM 

SlJMi=SM»'l+PCPEF<  I  ) 
3500    CONTINUE 
C  MHVI    cnR    SIMPLE    PEAL    COEcFICIENTS 

IF( IMAG.ME.2)    GO    TO    3504 

ZZDNM=ZDMM(1) 

ZZNUM=ZNUM( I) 

ZRDNM=ZPMM( l)*Zn( 1) 

QPFF(D=XNU^(1  )/X0NM(l  ) 

DOFFc( 1)=° 'TA**M/XnNM< 1) 

PO-Fc(  l  )  =11^-5  ^TE*P.ETA**M/XRDNM(  1) 

E0EFF(1)=XNIJM(1  )  /XROMM(l) 

SUM1=SUM1+PZF«=(  1  ) 
r  PR=LIMINARIES    TO    L DOPING    COR    TIM= 

3504    mjMMY=0.0 

PMNm=PrTA/(*ET4+DFLOAT(  M)*uPRAT5) 

ACURCY=RMND1+SUM1 

WR  I  TE  (  6  , 1 ) 
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FnoMATCl"  flOX,fCH«=CK    c0p 
WR*TE(6,2)    ACUPC.y,tcr 
=nqMAT<i     • , 10X t • TNDTCATnc 
1    FLAGGING    VALIT 
W?TT<-(6,3) 
FORMATC     '  ,11X 


=  » 


13/////) 
MEN3X,  'AVAIL 


ACCURACY    GF    NUMERICAL    SOLUTION'///) 

cOP    ACCURACY=«,~16.  12,  /  /,  11  X,  •  SUBRGUTI »  - 


BILITY',  5X,  *\ VAIL  /DOWN'  ,6X,  «E(  RE kir 


rLxn;ir       Mi^'l   n;',  JA,  '  AVfllLASll  1  I  Y'f  DA,  MV^  1L  /U'JW^  ■ 

1WAL1  ■  t7X,«F(TI  ME    ON)'//) 

RCNST=MPRATF*BrTA/(  BcTi  +  OFLOAT(  M)*IJDRATE) 

PCQNST=UFRATF*O=LOAT(M)*(Bi=T4+0  .5D0*UPR  ATE  *(  OFLOAK  M  )-  1 
1(R5TA+DFLQAT(M)*IJPRAT*:)**2 

ECNST=MPPATE*DFLOAT(M)*(OFLOAT(V 1+1.000) /( 2. 0D0*( BHTA+U 
1AT(M) )**2) 


1AT(M) )**2) 

DO    4500    L=l,100 
T-TSTPP*CFLCAT (L) 

RMMO=0. oro 
T~PM=0. 0^6 
OPMND=0 .000 
OTEPM=O.0PO 

FPMfjnrrfi  .COO 

RTEPM=o.OOO 

EPMND=0.00  0 

ETcRM=0 .000 

0'-    5000    I=IMAG,ILIM,2 

ZZR=ZP ( T) 

Q0S=DC0S(XR(2)*T) 

SIN=OST  N(X"  121  *T1 

"  -OFyr( xp (  1 )*T) 


000) )/ 
UPRATE*DFLO 


5000    COMTINMF 

IF(IMAG.rH.2J 
ReXPON=nF 


5001 


5003 
4500 


XP( (ROOTR(l) 
TERM=QFF  =  (  1  |«P«=X«>ON 

DTFR»=noF==<i  )*=-XPON 
RTrrRM  =  '?OrF  =  (l)*"FXPnN 
FTc-M=  ror*='(  l  )  *RFXFON 
-+T~" 


5001 

3ETA)*T> 


R 
ET. 

pp-iqyp  =  o .v,\oi  + 


T  op  qijp-q  v»n  l  +  'jD  MMr)+Qypp,M 

o  NV1  Al '  =  -F  r  HN  ST+R  r.NST*T+RR  MND+P  TER  M 

•=  X  PHP  =  P  C .  K'S  T  +  P  M  xlnl  *T  +  '  p  MN  0+  FT  CRM 

WCITF( 6, 5C03)     T,PPnBUPtOPRBUP,PNWAL»EXPUP 

FORMAT ( •     «.10X,F6.3,4X,P8.6.6X,F8.6,7X,F9.6,9X,C8. 

CONTINUE 


CHNTTNU': 

ST°c> 

END 


6) 
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APPENDIX   F:       SAMPLE    RESULTS 
CHFC*    C0D    ACCURACY    OF    NUMERICAL    SOLUTION: 

INHTCiTOR    -OR    ACCURACY=      1.000000000000 
SUBROUTINE    FLAGGING    V*L'JE=  0 


TIMJ 


iVAILABILITY 


AVAIL /DOWN 


E(PENFWAL) 


F(TIMF  ON) 


0.010 

0.997780 

0.020 

0.995  565 

0.03  0 

0.9^3356 

0.040 

0.991151 

0.053 

0.^88951 

0.060 

0.986756 

0.070 

0.984567 

0.080 

0.982384 

C.090 

0.980207 

0.100 

0.978038 

0.110 

0. C75S77 

0.120 

0.973725 

0.130 

0.971584 

0.14." 

0. 96°455 

0.150 

0.967339 

C.  160 

0.Q65238 

0.170 

0.963154 

0.180 

0.961089 

0.190 

0. 959044 

0.20  0 

0.957021 

0.210 

n .  9^n?^ 

0.220 

C. 953050 

0.230 

0.951  106 

0.24  0 

C.S49192 

0.250 

0.94730° 

0.26  0 

C. 945^60 

0 .  27  0 

0. 943645 

0.280 

0.941867 

C.290 

C.  94012.7 

0.300 

C. 938425 

0.310 

0.936^64 

0.320 

C. 935144 

0.330 

0.933566 

0.^.4  0 

0. 932031 

0.350 

C  930540 

0.36O 

0  .  Q2°OQ2 

0.370 

C. 927683 

0.38° 

C. 926329 

o  .300 

0.925015 

0.400 

C. 923745 

0.410 

0.92252n 

0.420 

0.921338 

0.430 

C.9202C1 

C.440 

0  .9i°107 

0.450 

0. 918C56 

0-460 

0.917C47 

0.470 

0.916030 

0.480 

0.915153 

0.49  0 

0.914267 

0.500 

0.913420 

0.000000 
0.0000C2 
0.000016 
0.00CO61 
o. 000172 
Q. 000394 
0.O00783 
0.001407 
0.002335 
0.003645 
0.005410 
0.007706 
0.010603 
C. 014166 
0.018453 
0.023512 
0.C29386 
0.036106 
0.C43692 
0.052157 
0.061504 
0.071727 
0.032810 
0.094731 
0.10~7461 
0.120962 
0.135194 
0.  150110 
0.165660 
0.1 8170Q 
0.198444 
0.215565 
0.233093 
0.25<~Q69 
0.269135 
0.287531 
O.306IOO 
0.32478"' 
0.343536 
0.362297 
0.381019 
0.39Q655 
n. 418162 
C. 436497 
0.454622 
0.472  502 
0.490103 
0.5^7396 
0.524355 
C. 540956 


0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000001 
0.000002 
0.00000  5 
0. 000009 
0.000015 
0.000025 
0. 000040 
0.000060 
0.000087 
0.0G0123 
0.00O170 
0.000229 
0.000301 
0.000390 
0. 000496 
0.0 0  06 2 2 
0, OOC779 
0.000942 
0.001139 
0.001363 
0.001617 
0.001901 
0.002218 
0.002569 
0.002°55 
0.003377 
0.003837 
0.004336 
0.004873 
0.005451 
0.006070 
0.006729 
0.007430 
0.C08173 
0.008957 
0.009783 
0.010650 
0.011559 
0. 012509 
0.013499 
0.0  145  29 
0.015599 
0.016707 
0.017854 
0.019037 


0.009989 
0.019956 
0.029900 
0.039823 
0.049723 
0.059602 
0.069458 
0.079293 
0.089106 
0.098897 
C. 108667 
0. 118415 
0.128141 
0.1378^7 
0.147531 
0.  157193 
0.166835 
0. 176^57 
0.186057 
0.  195638 
0 . 20  5198 
0.21473  8 
C. 224259 
C. 233760 
0.243243 
0.252707 
0.262152 
0.271580 
0.  280QC>0 
6.290382 
0.299753 
0.309118 
0.318461 
0.327789 
0.337102 
0.346400 
0.355634 
0  .364954 
0.374211 
0.38345^ 
0.392686 
0.401905 
0.411113 
0. 420309 
0.42°495 
0.438670 
0.44783  6 
C 456992 
0.4661^° 
0.475278 
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